% 
% close all;
clear all;
%% PLEASE ONLY CHANGE THE VALUES OF CAPITALIZED VARIABLES
NOIZE_LEVEL = 0.1;
MODE = 100;
LATTICE_SIZE = 200;
DX = 0.2;
%% 
m = 1:MODE;
z_phi = 1./(m.^2).*(1+NOIZE_LEVEL*randn(size(m)));
z_psi = 1./(m.^2).*(1+NOIZE_LEVEL*randn(size(m)));
%z_phi(1:end)=0;
%z_psi(2:end)=0;
figure(1);
subplot(2,1,1);
plot(m,z_phi);
title('Spectrum');
xlabel('mode');
ylabel('z_{\phi,n}');
subplot(2,1,2);
plot(m,z_psi);
xlabel('mode');
ylabel('z_{\psi,n}');
%%
% Generating surface
x_max = DX*(LATTICE_SIZE-1);
x = 0:DX:x_max;
h = zeros(size(x));

figure(2);
for n = 1:MODE
%    phi = phi_n(x,n);
%    psi = psi_n(x,n);
    [phi,psi] = eignfun_n(x,n);
    h = h+phi.*z_phi(n)+psi.*z_psi(n);
end
plot(x,h);

%%
% surface = surface1D(x,h);

%% Calculate m^2 from mode
K = (4/(x_max*DX*DX)).*((sin(pi*m/LATTICE_SIZE)).^2);
figure(3);
plot(m,K);
xlabel('mode');
ylabel('K');
m2_mode = sum(K.*(z_phi.^2+z_psi.^2));
%% Calculate m^2 by definition
slope = [h(1)-h(end),diff(h)]./DX;
% m2_surface = var(slope,1);
m2_surface = sum(slope.^2)/length(slope);
%%
% Calculate 
fprintf(1,'m^2 based on mode   : %.10f\n',m2_mode);
fprintf(1,'m^2 based on surface: %.10f\n',m2_surface);

%% Calculate r^2 from mode
r2_mode = (1/x_max)*sum(z_phi.^2+z_psi.^2);

%% Calculate r^2 from mode
r2_surface = var(h,1);
fprintf(1,'r^2 based on mode   : %.10f\n',r2_mode);
fprintf(1,'r^2 based on surface: %.10f\n',r2_surface);
